GIS with R: examples from Slovakia

Geospatial analysis. Why are we doing this?

1. Most of our data we collect are a recording of some phenomena that happened at a specific place and time.
2. Good visualisations and methods from spatio-temporal analysis field can help us understand our data.

This tutorial

1. Spatial data and coordinate systems.
2. Introduction to sf package.
3. Visualisations.
4. Example of analysis: spatial autocorrelation.

Extra resources

Recommended textbooks:

• Geocomputation with R. Available online.
• Wikle, Christopher K., Andrew Zammit-Mangion, and Noel Cressie. 2019. Spatio-Temporal Statistics with R. 1st ed. Boca Raton, Florida : CRC Press, [2019]: Chapman and Hall/CRC. https://doi.org/10.1201/9781351769723.

Other sources:

• r-spatial.org is a portal for people interested in spatial modelling.
• Packages like sf, sp tend to be well-documented, and with worked examples.

Spatial data

• Data frames, but now each row has a geometry associated with it (usually a point or an area).
• We need to be careful with the coordinate systems:
  • If working at macro scale (at world level), we need to take into account curvature of the Earth.
  • We stick to coordinate systems where we can work with Euclidean geometry.
  • The cadaster office in Slovakia uses krovak system, which is the legally binding coordinate system (coded as EPSG:5513).
  • This website is useful for reading off the coordinates for a particular point.

Obtaining geometry (map) files

• The map files come in different formats. One of the most common ones are shapefiles (.shp), geodatabase (.gdb), geopackage (.gpk), geoJSON (.geojson).
• This page in the geoportal.sk has the necessary files in different formats. We stick to the shapefiles (.shp) format. The portal allows downloading the maps at different resolutions. In this guide, we will be working with the most precise one, which is referred to as ‘Základná úroveň/ ZBGIS - Administratívne hranice’.

sf package

Install the package using:

install.packages('sf')

library('sf')

st_read() to load maps

st_read(dsn, layer)

• dsn: file name of the map file.
• layer: certain file formats allow saving multiple layers, we will work with single layer files (it does not mean we can’t have multiple indicators/statistic for each location).

Load the map of Slovakia at the municipality level:

mapa.obec.sf <- st_read(dsn = file.path(getwd(),
                                        'maps', 
                                        'ah_shp_0', 
                                        'obec_0.shp'))

Inspecting the map object - coordinates

To see what coordinate system it use st_crs() function:

st_crs(mapa.obec.sf)

## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813975277778 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs"

Due to the missing code, we do a transformation to the official code using st_transform() function.

crsCodeSkCadastreOffice <- 5513  # save the cadastre CRS to a constant
mapa.obec.sf <- st_transform(x = mapa.obec.sf, 
                             crs = crsCodeSkCadastreOffice)
st_crs(mapa.obec.sf)                   # 

Inspecting the map object - data

summary(mapa.obec.sf)

##       DOW                  AUT         ACH         SOI       FACC           IDN4              NM4            IDN3      
##  Min.   :2019-04-30   Min.   :1   Min.   :1   Min.   :9   FA004:2927   Min.   :500011   Lúčka   :   4   Min.   :101.0  
##  1st Qu.:2019-04-30   1st Qu.:1   1st Qu.:1   1st Qu.:9                1st Qu.:508989   Porúbka :   4   1st Qu.:402.0  
##  Median :2019-04-30   Median :1   Median :1   Median :9                Median :518204   Dúbrava :   3   Median :606.0  
##  Mean   :2019-04-30   Mean   :1   Mean   :1   Mean   :9                Mean   :521206   Lúčky   :   3   Mean   :547.7  
##  3rd Qu.:2019-04-30   3rd Qu.:1   3rd Qu.:1   3rd Qu.:9                3rd Qu.:526138   Petrovce:   3   3rd Qu.:708.0  
##  Max.   :2019-04-30   Max.   :1   Max.   :1   Max.   :9                Max.   :599972   Píla    :   3   Max.   :811.0  
##                                                                                         (Other) :2907                  
##               NM3            IDN2                    NM2          VYMERA            Shape_Leng       Shape_Area       
##  Košice - okolie: 114   Min.   :1.000   Prešovský      :665   Min.   :   361255   Min.   :  3211   Min.   :   354208  
##  Rimavská Sobota: 107   1st Qu.:4.000   Banskobystrický:516   1st Qu.:  7207245   1st Qu.: 14031   1st Qu.:  7207035  
##  Prešov         :  91   Median :6.000   Košický        :461   Median : 11655767   Median : 18582   Median : 11647728  
##  Levice         :  89   Mean   :5.415   Nitriansky     :354   Mean   : 16752332   Mean   : 21367   Mean   : 16747701  
##  Bardejov       :  86   3rd Qu.:7.000   Žilinský       :315   3rd Qu.: 19769322   3rd Qu.: 25503   3rd Qu.: 19794868  
##  Trebišov       :  82   Max.   :8.000   Trenčiansky    :276   Max.   :359784254   Max.   :204284   Max.   :359936992  
##  (Other)        :2358                   (Other)        :340                                                           
##           geometry   
##  MULTIPOLYGON :2927  
##  epsg:5513    :   0  
##  +proj=krov...:   0  
##                      
##                      
##                      
## 

Inspecting the map object - data

From this, it is obvious that:

1. IDN4 is a unique identifier for each municipality.
2. NM4 is the name of the municipality.
3. IDN3 is the identifier for okres.
4. NM3 is the name for okres.
5. IDN2 is the identifier for kraj.
6. NM2 is the name for kraj.
7. Shape_Leng and Shape_Area are the perimeter and the area of the municipality.
8. It’s not clear to me what VYMERA is given that the values slightly deviate from the Shape_Area.

The other columns are not relevant for this tutorial.

In the visualisation section, we will stick to these values even though they are not particularly interesting. In later sections we will add more interesting data.

Data preparation

Before, we proceed, we rename the columns to more meaningful names:

library(dplyr)
mapa.obec.sf <- mapa.obec.sf %>%
  rename(obec = NM4) %>%
  rename(obec_id = IDN4) %>%
  rename(okres = NM3) %>%
  rename(okres_id = IDN3) %>%
  rename(kraj = NM2) %>%
  rename(kraj_id = IDN2) %>%
  rename(obvod_uzemia = Shape_Leng) %>%
  rename(plocha_uzemia = Shape_Area)

Visualisation - options

Some of the most popular options are:

1. default plot() method in the sf package,
2. ggplot,
3. tmap.

We will give examples of all throughout the rest of this tutorial.

Example: migration in Slovakia I

1. Load the migration and population data.
2. Merge the data with the map files.
3. Visualise with migration type breakdown:
   • within-country migration
   • between-countries migration

Example: migration in Slovakia II

Load migration and population data

migration.district.df <- read_excel(file.path(getwd(), 'data', "raw", "migration_by_district.xlsx"),
                                    skip=8, col_names=FALSE)
colnames(migration.district.df) <- c('district', 'inflow_sk', 'inflow_abroad', 'outflow_sk', 'outflow_abroad')
migration.district.df <- migration.district.df %>%
  dplyr::mutate(district = gsub("^District of (.+)$", "\\1", district))

# load population data and clean them
population.municipality.df <- read_excel(file.path(getwd(), 'data', "raw", "population_municipality.xlsx"),
                                         skip=7, col_names=FALSE)
colnames(population.municipality.df) <- c('district', 'municipality', 'population')
population.district.df <- population.municipality.df  %>%
  dplyr::mutate(district = gsub("^District of (.+)$", "\\1", district)) %>%
  dplyr::group_by(district) %>%
  dplyr::summarise(population = sum(population))

Example: migration in Slovakia II

Merge the data with maps

library(stringi) # for removing diacritis for joins

pop.migration.district.df <- migration.district.df %>%
  dplyr::full_join(population.district.df, by = 'district') %>%
  dplyr::mutate(district_ascii = stri_trans_general(str = district, id = "Latin-ASCII"))

pop.migration.district.sf <- mapa.okres.sf %>%
  dplyr::mutate(district_ascii = stri_trans_general(str = okres, id = "Latin-ASCII")) %>%
  dplyr::left_join(pop.migration.district.df, by = "district_ascii")

pop.migration.district.sf <- pop.migration.district.sf %>%
  dplyr::mutate(outflow_abroad_percent = 100*(outflow_abroad / population)) %>%
  dplyr::mutate(outflow_sk_percent = 100*(outflow_sk / population)) %>%
  dplyr::mutate(inflow_abroad_percent = 100*(inflow_abroad / population)) %>%
  dplyr::mutate(inflow_sk_percent = 100*(inflow_sk / population))

Plot using the default plot() method

plot(x = pop.migration.district.sf[c('outflow_abroad_percent')], 
     key.pos = 1, main='% of population moved abroad')

Example: migration in Slovakia III

Example: migration in Slovakia IV

Plotting side-by-side using the default method is possible, but has limitations such as not supporting legends. tmap package can handle those very well.

This library allows integration with Open Street Map data and interactive maps on top of regular plotting. What we show here is only scratching the surface, for more examples, please visit the tmap’s github page.

install.packages('tmap')

library(tmap)

Example: migration in Slovakia V

pop.migration.sf.tall <- pop.migration.district.sf %>%
  tidyr::gather("type", "percent", c("outflow_abroad_percent", "outflow_sk_percent",
                                     "inflow_abroad_percent", "inflow_sk_percent")) %>%
  dplyr::mutate(type = gsub("_", " ", gsub("(.+)\\_percent", "\\1", type)))

sk.migration <- pop.migration.sf.tall %>%
  dplyr::filter(grepl('.*sk', type))

tm_shape(sk.migration) +
  tm_fill(col = "percent", style = "cont", title="% of population") + 
  tm_style('gray', title = 'Migration within Slovakia') + 
  tm_facets('type')

Exercise: do the same for migration to/out of Slovakia.

Example: migration in Slovakia V

Example: deaths in Slovakia

here show ggplot

• Download the dataset from the official statistics portal
• Do some plots, by sex…

Example: death causes in Slovakia

• Download the dataset from the official statistics portal
• Break down by disease categories

Using spatial analysis in regression models

The example we give is very crude, but we just want to illustrate the point.

• check if there is a spatial correlation for the rates of deaths / age at which people die
• fit some glm model and see if the residuals are correlated…

here show tmap methods

1. Check spatial correlation for death rate for plucne choroby
2. Try to fit the thing about altitude stuff….
3. Fit number of breathing deaths a function of altitude
4. check residuals…

Visualisation - default plot() method

plot(x = mapa.okres.sf[c('plocha_uzemia')],
     key.pos = 1) # this adds the legend to the bottom (1)

Visualisation - default plot() method

This is a good place to start when doing exploratory analysis, but it’s better to use more sophisticated functions to produce better-looking plots.

Visualisation - ggplot

ggplot can also be used in conjunction with spatial objects.

library(ggplot2)  # Load the ggplot library
library(viridis)  # Load the library for a nice colour palette
ggplot(mapa.okres.sf) +
  geom_sf(aes(fill = plocha_uzemia)) +
  scale_fill_viridis("Plocha m\u00B2") +
  ggtitle("Plocha okresov") +
  theme_bw()

Visualisation - ggplot

Visualisation - tmap

This is a very versatile library that allows integration with Open Street Map data and interactive maps on top of regular plotting. What we show here is only scratching the surface, for more examples, please visit the tmap’s github page.

Before we proceed, we install the libraries and load them.

install.packages('tmap')

library(tmap)

Visualisation - tmap

## tmap mode set to plotting

tmap_mode("plot")
tm_shape(mapa.okres.sf) +
  tm_fill(col = "plocha_uzemia", style = "cont", title="Plocha m\u00B2") + 
  tm_style('gray', title = 'Plocha okresov m\u00B2')

Visualisation - tmap

Visualisation - tmap

For example, we can do a breakdown by kraj.

tm_shape(mapa.okres.sf) +
  tm_fill(col = "plocha_uzemia", style = "cont", title="Plocha m\u00B2") + 
  tm_style('gray', title = 'Okresy SR') + 
  tm_facets('kraj')

Visualisation - tmap

Visualisation - tmap

For plots with fewer components (e.g. kraj), we might want to add labels. Be careful with the overlaps of the labels and overflows outside the frame. This post gives a good solution for the latter.

tm_shape(mapa.kraj.sf) +
  tm_fill(col = "plocha_uzemia", style = "cont", title="Plocha m\u00B2") + 
  tm_style('gray', title = 'Okresy SR') + 
  tm_text('kraj', auto.placement = TRUE)

To save the map, use tmap_save() function.

Visualisation - tmap

mapa.okres.sf.gps <- st_transform(x = mapa.okres.sf, crs=4326)
osm_tiles = tmaptools::read_osm(st_bbox(mapa.okres.sf.gps), type = "osm")
tmap.object <- qtm(osm_tiles, raster.alpha = 0.6) + 
  tm_shape(mapa.okres.sf.gps) +
  tm_fill(col = "Shape_Area", style = "cont", title="Plocha m\u00B2") + 
  tm_style('gray', title = 'Okresy SR')

A few useful visualisations

Case study: presidential elections 2019

1. Download the raw dataset of the results from here.
2. Load the results, merge it with the geography you are interested in (obec, kraj)
3. Plot the results.

Case study: presidential elections 2019

Load the data:

library(readxl)
results.obec.file.path <- file.path(getwd(), 'data', "raw", 
                                    "PRE_2019_KOLO1_xlsx",
                                    "PRE_2019_KOLO1_tab03e.xlsx")
results.obec <- read_excel(results.obec.file.path)
results.obec <- results.obec[-(1:4), ]
names(results.obec) <- c('kraj_id', 'kraj', 'uzemny_obvod_id', 'uzemny_obvod',
                         'okres_id', 'okres', 'obec_id', 'obec',
                         'candidate_id', 'candidate_fname', 'candidate_sname',
                         'platne_hlasy_count', 'platne_hlasy_percent',
                         'pullouts')
results.obec <- results.obec %>%
  dplyr::mutate_at(c("platne_hlasy_count", "platne_hlasy_percent"), as.numeric)

Case study: presidential elections 2019

Join the data with the map:

results.obec.sf <- mapa.obec.sf %>%
  dplyr::select(obec_id) %>%
  dplyr::mutate(obec_id = as.character(obec_id)) %>%
  dplyr::full_join(results.obec, by = c("obec_id" = "obec_id")) %>%
  dplyr::mutate(name = paste(candidate_fname, candidate_sname, sep = " "))

Case study: presidential elections 2019

Since we are working with ‘obec’ level, we focus only on one ‘kraj’ region, and we only compare top 4 candidates in that region.

zilinsky.kraj.results.sf <- results.obec.sf %>%
  dplyr::filter(kraj_id == 5)  # zilinsky kraj

# select top 4 candidates in Zilinsky kraj
zilinsky.kraj.top.candidates  <- zilinsky.kraj.results.sf %>%
  dplyr::group_by(name) %>%
  dplyr::summarise(vote_total = sum(platne_hlasy_count)) %>%
  dplyr::arrange(vote_total) %>%
  dplyr::top_n(4, vote_total) %>%
  dplyr::pull(name)

# filter out only the 4 candidates
zilinsky.kraj.results.top.candidates.sf <- zilinsky.kraj.results.sf %>%
  dplyr::filter(name %in% zilinsky.kraj.top.candidates)

Case study: presidential elections 2019

tmap_mode("plot")
tm_shape(zilinsky.kraj.results.top.candidates.sf) +
  tm_fill(col = "platne_hlasy_percent", style = "cont", title="Podiel hlasov %") + 
  tm_style('gray', title = 'Žilinsky kraj') + 
  tm_facets('name')

Case study: presidential elections 2019

Case study: presidential elections 2019

As we can see this visualisation is not ideal, so we go for the interactive map where we can zoom in/out as we please.

library(tidyr)

zilinsky.kraj.results.top.2.sf <- zilinsky.kraj.results.top.candidates.sf %>%
  dplyr::filter(name %in% c('Maroš Šefčovič', 'Zuzana Čaputová'))

zilinsky.results.wide.sf <- zilinsky.kraj.results.top.2.sf %>%
  dplyr::select(name, platne_hlasy_percent, obec) %>%
  tidyr::spread(name, platne_hlasy_percent)

tmap_mode("view")
tm <- tm_shape(zilinsky.results.wide.sf) + 
  tm_polygons(c('Maroš Šefčovič', 'Zuzana Čaputová')) + 
  tm_facets(sync = TRUE, ncol = 2) +
  tm_layout(legend.text.size = 0.4)
tm

Case study: presidential elections 2019